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Abstract 

In a series of papers, Olsson (1994, 1995), Olsson & Oliger (1994), Strand (1994), Gerritsen 
& Olsson (1996), Yee et al. ( 1999a, b, 2000) and Sandham & Yee (2000), the issue of nonlinear 
stability of the compressible Euler and Navier-Stokes Equations, including physical boundaries, 
and the corresponding development of the discrete analogue of nonlinear stable high order 
schemes, including boundary schemes, were developed, extended and evaluated for various fluid 
flows. High order here refers to spatial schemes that are essentially fourth-order or higher away 
from shock and shear regions. The objective of this paper is to give an overview of the progress 
of the low dissipative high order shock-capturing schemes proposed by Yee et al. ( 1999a, b, 
2000). This class of schemes consists of simple non-dissipative high order compact or non- 
compact central spatial differencings and adaptive nonlinear numerical dissipation operators 
to minimize the use of numerical dissipation. The amount of numerical dissipation is further 
minimized by applying the scheme to the entropy splitting form of the inviscid flux derivatives, 
and by rewriting the viscous terms to minimize odd-even decoupling before the application of 
the central scheme (Sandham & Yee). 

The efficiency and accuracy of these schemes are compared with spectral, TVD and fifth- 
order WENO schemes. A new approach of Sjogreen & Yee (2000) utilizing non-orthogonal 
multi- resolution wavelet basis functions as sensors to dynamically determine the appropriate 
amount of numerical dissipation to be added to the non-dissipative high order spatial scheme 
at each grid point will be discussed. Numerical experiments of long time integration of smooth 
flows, shock-turbulence interactions, direct numerical simulations of a 3-D compressible tur- 
bulent plane channel flow, and various mixing layer problems indicate that these schemes are 
especially suitable for practical complex problems in nonlinear aeroacoustics, rotorcraft dy- 
namics, direct numerical simulation or large eddy simulation of compressible turbulent flows 
at various speeds including high-speed shock-turbulence interactions, and general long time 
wave propagation problems. These schemes, including entropy splitting, have also been ex- 
tended to freestream preserving schemes on curvilinear moving grids for a thermally perfect 
gas (Vinokur & Yee 2000). 


I. Motivation and Overview 

Strong theoretically based high order high-resolution shock-capturing schemes have domi- 
nated algorithm development for fluid flows for nearly tw'o decades. During the 1990’s, the 
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focus was on increasing the order of accuracy for the interior grid points of the computa- 
tional domain with stable numerical boundary conditions generally not included as part of the 
theoretical development. Multidimensional upwinding and high order unstructured or finite 
element constructions have flooded the computational fluid dynamics (CFD) literature. The 
majority of these schemes are either too CPU intensive for practical applications or they 3se 
still in the development stage. In spite of their high-resolution capability for rapidly evolving 
flows and short time integrations, these schemes often exhibit undesirable amplitude errors 
for long time integrations. On the other hand, simplicity, efficiency and highly parallelizable 
robust algorithms are the major requirements in industrial, geophysical, space exploration and 
military practical CFD applications. The objective of this paper is to give an overview of the 
progress of a class of simple, highly parallelizable high order shock-capturing schemes that 
meets many of the requirements for practical computations, especially for long time integra- 
tions of unsteady flows. For ease of reference, “schemes” or more precisely “interior schemes” 
here generally refer to spatial difference schemes for the interior grid points of the computa- 
tional domain, whereas “boundary schemes” are the numerical boundary difference operators 
for grid points near the boundaries. However, without loss of generality, we also adopt the 
conventional terminology of denoting “scheme” as either the “combined interior and boundary 
schemes” or just the "interior scheme” interchangeably within the context of the discussion. 

Before 1994, rigorous stability estimates for accurate and appropriate boundary schemes as- 
sociated with fourth-order or higher spatial interior schemes were the major stumbling block in 
the theoretical development of combined interior and boundary schemes for nonlinear systems 
of conservation laws. Spatial nonlinear stability of initial boundary value problems (IBVPs) 
goes hand-in-hand with the appropriate amount of nonlinear numerical dissipation required to 
stabilize the combined interior and boundary schemes. The delicate balance of the numerical 
dissipation for stability without the expense of excessive smearing of the flow physics after long 
time integrations poses a severe challenge for unsteady flow simulations of this type. Actually, 
there are two possible sources of stabilizing mechanisms; namely, (a) from the governing equa- 
tion level and (b) from the numerical scheme level. Employing a nonlinear stable form of the 
governing equations 'more conditioned form of the PDE) in conjunction with the appropriate 
nonlinear stable scheme for IBVPs is one way of minimizing the use of numerical dissipation. 

The major tool used to overcome the stumbling block is a generalized energy method. The 
basic building block in establishing a stable energy estimate for high order spatial central 
schemes for nonlinear hyperbolic conservation laws relies on the aforementioned mechanisms 
(a) and (b). From the governing equation level, a special transformation of the conservation 
laws to an appropriate form for the application of the continuous energy estimate for a stable 
IBVP of the governing equations is needed. From the numerical scheme level, a compatible 
boundary scheme for high orrler central interior schemes that satisfy the discrete analogue of 
the continuous energy' estimate is needed. See Strand (1994), Olsson (1995), Olsson fe Oliger 
(1994) and references cited. 

Olsson proved that an energy estimate can be established for second-order central schemes. 
To obtain a rigorous energy estimate for high order central schemes, one must apply the 
scheme to the split form of the inviscid governing equation. For the Euler equations, the 
transformation consists of a convex entropy function that satisfies a mathematical entropy 
condition. This mathematical entropy function, in this case, can be a function of the physical 
entropy. Therefore, the resulting splitting is hereafter referred to cis entropy splitting for ease 
of reference. Here, the entropy splitting should not be confused with the traditional flux vector 
splittings such as the Steger and Warming splitting (1981) or other variants. The traditional 
flux vector splitting splits the flux function into different parts and most often into upwind 
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and downwind portions. However, the entropy splitting splits the flux derivatives and the 
time derivative of the conservative dependent variables using the properties of the chosen 
mathematical entropy function and the symmetrizability of the conservation laws, without 
reference to any upwinding. For the viscous terms, only symmetry is needed in the derivation 
of the energy estimate. These developments have made high order non-dissipative spati^d 
central schemes of renewed interest to CFD practitioners where simplicity, efficiency and high 
parallelizability are their treidemark. 

The recently developed high order low-dissipation shock capturing schemes using character- 
istic filters of Yee et al. (1999a) fit in the entropy splitting framework. Instead of applying a 
scalar dissipation or filter (Gerritsen & Olsson), they supply nonlinear filters based on, locally, 
the different wave characteristics of the inviscid flux. For complex shock waves, shear and 
turbulence interactions, one has better control of the amount of dissipation associated with 
each wave. For efficiency, Yee et al. proposed a combination of a high order base scheme and 
a nonlinear filter operator. The base scheme consists of narrow grid stencil high order compact 
or non-compact centered non-dissipative classical spatial differencings. The filter consists of a 
product of the dissipative portion of a low order total variation diminishing (TVD), essen- 
tially non-oscillatory (ENO) or weighted ENO (WENO) scheme and an artificial compression 
method (ACM) sensor. In contrast to hybrid schemes that switch between spectral or spectral- 
like non-shock-capturing schemes and high order ENO schemes, the high order non-dissipative 
base scheme is alw^ays activated. The role of the ACM sensor is to reduce the amount of nu- 
merical dissipation away from shock and shear regions. As an alternative to the ACM sensor, 
Sjogreen and Yee (2000) utilized non-orthogonal wavelet basis functions as multi-resolution 
sensors to dynamically determine the amount of nonlinear numerical dissipation to be added 
at each grid point. The resulting sensor function is also readily usable for grid adaptation 
purposes. The final grid stencil of these schemes is five points in each spatial direction if 
second-order TVD schemes are used as filters, and seven points if second-order ENO schemes 
are used as filters for a fourth-order base scheme. Studies showed that higher accuracy was 
achieved with less CPU time and fewer grid points when compared with that of standard high 
order TVD. positive. ENO or WENO schemes. Table 1 shows the flow chart of the schemes. 

The studies in Yee et al. (1999b, 2000) and Sandham & Yee indicate that entropy splitting 
can improve the overall stability of the scheme, and that the amount of numerical dissipation, 
if needed, is less than for the unsplit approach. They view entropy splitting as a conditioned 
form of the original governing equations. Here, “condition the governing equation” is different 
from “preconditioning of the flow equations or their discretized counterparts” in convergence 
acceleration of time- marching to steady states. Their studies also indicate that entropy splitting 
alone can improve nonlinear stability even when one employs boundary schemes that do not 
satisfy the discrete generalized energy estimate. This stability property of the entropy splitting 
is valuable not just for the class of schemes in question, but can also be applied to other schemes 
commonly used in practical CFD applications. This emphasizes the fact that one should always 
try to apply numerical schemes to a more conditioned form of the governing equations. 

Extension of these schemes to freestream preserving schemes for 3-D curvilinear moving grids 
for a thermally perfect gas is reported in Vinokur and Yee (2000). The main difficulty in the 
extension of high order schemes to curvilinear grids is the high order numerical evaluation of the 
geometric terms, arising from the coordinate transformation, to satisfy a coordinate-invariant 
freestream preservation condition. The question of the extendibility of the entropy splitting 
concept to other physical ecjuations of state and evolutionary equation sets was examined 
in Yee et al. (1999b. 2000). Their study shows that the entropy splitting can be formally 
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extended to a thermally perfect gas, with the internal energy being an arbitrary function 
of temperature. For non-equilibrium flows which consist of a mixture of different species, 
each obeying a thermally perfect gas law, extension of the splitting is problematic. While 
they were able to prove the symmetry and homogeneity properties required for the energy 
estimate, the degree of homogeneity can only be obtained by solving a system of nonlinear 
equations. In addition, to obtain the Jacobian of the transformation required inverting a 
non-sparse linear system. It would therefore be difficult to establish the positive definite 
condition in closed form. Consequently, the extension of the method to non-equilibrium flows 
is not practically feasible. If the homogeneity condition is not required, then one can use 
symmetry variables based on the physical entropy, as was shown by Chalot et al. (1990). In 
this case, the resulting PDFs are in pure non-conservative form and entropy splitting is no 
longer applicable. For magnetohydrodynamics (MHD), the magnetic field has to be added 
as a ‘conservative'’ variable. But the square of the magnetic field is one of the terms in the 
definition of the total energy-. Thus, from dimensional arguments it is clear that one cannot 
obtain the homogeneity condition. A similar situation exists for the artificial compressibility 
method of treating incompressible flow. In the Maxwell equations, we have a linear system 
of hyperbolic equations that can easily be symmetrized. Thus Strand’s numerical boundary 
operators are still valid, but entropy splitting is not needed. For non-equilibrium flows, if one 
solves the species and flow equations separately in a loosely-coupled manner, then the flow 
equations effectively satisfy a locally thermally perfect gets law and a '‘local” form of entropy 
splitting is applicable. 

We would like to point our that although the formal extension of entropy splitting is lim- 
ited to a thermally perfect gas, the numerical schemes themselves do not have this restric- 
tion. Consequently, the schemes discussed here are applicable to equilibrium real gas, non- 
equilibrium and the artificial compressibility method of treating incompressible flows, MHD 
and the Maxwell equations. 

II. Entropy Splitting and Numerical Methods 

This section reviews the entropy splitting formula for the 2-D compressible Euler equations 
for a perfect gas in Cartesian coordinates. Formulas for the corresponding 3-D case can be 
found in Appendix B of Yee et al. (1999b) and for curvilinear moving grids for a thermally 
perfect gas in Vinokur and Yee (2000). The mathematical theory is quite involved. Interested 
readers are referred to references cited. The Yee et al. (1999a) and Sjogreen & Yee numerical 
methods used in conjunction with the entropy splitting are also summarized. 

2.1. Summary of Entropy Splitting for a Perfect Gas 

In vector notation the 2-D compressible time-dependent Euler equations in conservation 
form for an equilibrium real gas can be written, in Cartesian coordinates, as 

f/f + + = 0, (la) 

where Ut = Fj^ = ^ and Gy — ^ and the J7, F, G, are vectors given by 


■ p ' 


pu 


pv 

pu 

: F = 

pu^ + p 

: G = 

puv 
2 1 

pi- 


puv 


-f p 

. e . 


. eu 4- pu . 


.ev pv. 


(lb) 


5 


The dependent variable U is the vector of conservative variables, and {p,u,v,p) is the vector 
of primitive variables. Here p is the density, u and v are the velocity components, pu and 
pv cire the x- and y-components of the momentum per unit volume, p is the pressure, e = 
p[e + (u^ + u^)/2] is the total energy per unit volume, and £ is the specific internal energy. For 
a thermally perfect gas, the equation of state is p = pRT, where R is the specific gas constant, 
and T is the temperature with e = s{T). For constant specific heats (calorically perfect gas) 
£ = c„T, where c„ is the specific heat at constant volume. 

The eigenvalues associated with the flux Ja<;obian matrices of F and G are (u, u, u ± c) 
and {v,v,v ± c). where c is the sound speed. The two u, u and v, v characteristics are line 2 u:ly 
degenerate. Hereafter, we refer to the fields associated with the u±c and v±c characteristics 
£is the nonlinear fields and the fields associated with the u, u and v, v characteristics cis the 

linear fields. 

The entropy splitting for the compressible Euler equations for a perfect gas utilizes the 
result of Harten (1983) on symmetric form for systems of conservation laws as the backbone. 
The idea is to introduce a symmetry transformation from the vector of conservative variables 

to a new vector \V, referred to as the “entropy variables”. The transformation is chosen 
so that Fw = ^ and Gw = -§w symmetric, and Uw = ^ symmetric and positive 
definite. One then restricts the transformations to those that allow a special splitting of Fj^.Gy 
and Ut- This requires that the entropy variable W be chosen such that F{JJ{W)), G{U{W)) 
and U{W) are homogeneous functions of W of degree ,5. Then the splitting of results in 
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The corresponding U’ can be written as 


ir = — [e - 2pc - p(l + 5) -pu -pu p]^; P* - -{pp ^) • (3) 

P 

A similar splitting can be written for Gy and Ut- The forms for Fw and Gw can be found in 
Yee et al. (1999b). Under the above conditions one can rigorously establish a bound on the 
rate of growth of the energy norm in terms of the eigenvalues corresponding to the incoming 
characteristic variables at the boundary of the domain. 

Normally, we need to compute Uw for the split form of Ut = -j^Ut + -^^UwWt- However, 
we only consider a semi-discrete approach of applying temporal discretizations. Aside from 
using the split form of the inviscid flux derivatives Fj, and Gy. we do not have to use the 
split form of L't for implementation. Thus the final form of the semi-discrete entropy splitting 
approach still can be expressed in terms of conservative and primitive variables, making possible 
easy and efficient implementation in existing computer codes. The splitting parameter 5 has 
to satisfy 5 > 0 or 5 < Harten only considered 5 < This choice of 5 appears to be 
“non-standard" or "nonphysical” in the sense that more than 100% of the conservative portion 
and a negative non-conservative portion is used. Although Gerritsen and Olsson considered 
the 5 > 0 range which Harten overlooked, they set 5 = 1 >n conjunction with high order 
central differencing schemes in all of their numerical examples. This choice of 5 corresponds to 
splitting of the flux derivative into equal conservative and non-conservative proportions. Yee 
et al. (1999b. 2000) recommend the use of /? > 0. The degree of improvement in stability over 
the unsplit approach depends on the choice of 5- For 5 > 100 the benefit from the entropy 
splitting is diminishing, since this case is close to the unsplit situation. In addition, the choice 
of 5 is also problem- dependent. For certain problems, e.g., complex shock-shear interactions, 
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the beneficial range of j3 can be small. A combination of wavelets cis filter (Sjogreen & Yee) 
sensors and grid adaptations (Gerritsen & Olsson) might be able to enlarge the beneficial 
range since the resulting wavelet sensor function is also readily applicable for grid adaptation 
purposes. The grid adaptation might be able to minimize the underresolved grid related 
spurious oscillation producing instability. This will be a subject of future research, 

2.3. Numerical Methods 

The spatial discretizations considered in Yee et al. (1999a) consist of two paits, namely, 
a base scheme and a filter. When numerical dissipations or filters are not used, the scheme 
consists of only the base scheme. If entropy splitting is used, the base scheme is applied to the 
split form of the inviscid flux derivatives. Possible non-dissipative high order base schemes for 
Fjc and Gy and the viscous terms (if present) are the standard fourth and sixth-order compact 
and non-compact central schemes for the interior grid points. 

There are many possible candidates for the filter operator in conjunction with high order 
base schemes. For efficiency and ease of numerical boundary treatment, Yee et al. (1999a) 
proposed using filter operators whose grid stencils have a width similar to that of the base 
scheme. The filter operator consists of the product of a sensor and a nonlinear dissipation. See 
Table 1 for the road map. Denote Fj^k as the discrete approximation of the inviscid flux F at 
ijAx.kAy). where Ax and Ay are the grid spacings in the x- and y-directions and j and k 
are the corresponding spatial indices. Let the filter vector in the x-direction be of the form 






(^) 


F*^ 1 k modified form of the nonlinear dissipation portion of the standard numerical 

flux. For characteristic based methods, the quantity (wdth the k index suppressed) is the 

right eigenvector matrix of using Roe’s average state (Roe's approximate Riemann solver). 
We define G" , . in the same manner. The elements of , i (with the k index suppressed), 

j.A: + ^ J ' I 


denoted by (6^ are 

J-r 2 



( 5 ) 


in (5) is the dissipative portion of the high resolution scheme resulting from using a TVD, 
MUSCL with slope limiters. ENO or WENO scheme. Formulae for ,i are well known and 
can be found in the literature. See Yee et al. (1999a) for details and for a discussion of other 
possible filters. 

Here is the sensor and is a mechanism for controlling excess numerical dissipation that 

is inherent in the dissipative portion of standard high-resolution shock-capturing schemes. Two 
possible sensors are considered. They are the ACM sensor and the wavelet sensor (Sjogreen & 
Yee). 

ACM Sensor: For the ACM filter, == parameter k is problem-dependent. 

For smooth flows that are absent of high shears, k can be very small. It is used to minimize 
spurious high frequency oscillation producing nonlinear instability associated with pure central 
schemes, especially for long time integration problems. Different physical problems require 
different values of k because of the large variation in flow properties. The k value may vary 
from one characteristic wave to another, and from one region of the flow field to another region 
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with different flow structure. The function 0^ ^ is the Harten ACM gradient sensor but utilized 
in a different context than Harten originally intended. For a general 2m + 1 point base scheme, 
Harten recommended 



= max 
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Here the parameter p is an exponent > 1 and is not the “pressure p’’ in (1). Instead of varying 
K for the particular flow problem, one can vary p. For larger p, less numerical dissipation is 
added. Note that by varying p > 1 in (7), one can essentially increase the order of accuracy of 
the filter. The are elements of + Uj^k)- The corresponding 

using the MUSCL formulation are instead functions of the left and right states of {/. 

Multi- Resolution Non-orthogonal Wavelet Sensors: To avoid the tuning of the arbitrary 
parameters k and/or p in (7), one can replace i by a wavelet sensor , i (Sjogreen & 

Yee, 2000). With a proper choice of the wavelet basis function, we have a better control 
on the proper distribution of numerical dissipations leading to a more accurate simulation 
than the ACM sensor. Wavelets were originally developed for feature extraction in image 
processing and for data compression. It is well known that the regularity of a function can be 
determined from its wavelet coefficients (Daubechies 1992, Mallat & Zhong 1992) far better 
than from its Fourier coefficients. By computing wavelet coefficients (of a suitable wavelet 
basis function), we obtain very precise information about the regularity of the function in 
question. As of the 1990's, wavelets are a new' class of basis functions that are finding use in 
analyzing and interpreting turbulence data from experiments. They also are used for analyzing 
the structure of turbulence from numerical data obtained from DNS or large eddy simulation 
(LES). See Farge (1992) and her later work, and Perrier et al. (1999). Recently, wavelets have 
been used for grid adaptation (Gerritsen & Olsson) and to replace existing basis functions 
in constructing more accurate finite element methods. Here we utilize wavelets to adaptively 
control the amount of numerical dissipation that is inherent in standard high-resolution shock- 
capturing schemes. The resulting wavelet sensors are readily available as more desirable grid 
adaptation indicators than the commonly used grid adaptation indicators. 

The wavelet sensor of Sjogreen & Yee is obtained by computing the so called “Lipschitz 
Exponent” of a chosen vector to be sensed with a suitable multi-resolution non-orthogonal 
wavelet basis function that is capable of detecting shocks, shears, spurious oscillations and 
turbulence. Here, "vectors or variables to be sensed” means the represented vectors or variables 
that are suitable for the extraction of the desired flow physics. The study in Sjogreen &: Yee 
showed that for a proper choice of the wavelet basis function, the wavelet sensor is physical 
problem-independent for all of the test cases considered in Yee et al. ( 1999a, b). The variables 
to be sensed can be the density and/or pressure, the characteristic variables, the or 

the entropy variables W. There are two types of non-orthogonal wavelet basis functions that 
Sjogreen & Yee considered. One is similar to the B-spline wavelet (Mallat & Zhong) used by 
Gerritsen & Olsson for grid adaptation and the other is modification of the multi-resolution 
method of Harten (1995) as a redundant multi-resolution wavelet. The B-spline wavelet sensor 
requires slightly more arithmetic operations than the redundant form of Harten wavelet sensor. 
The fined form for involves mainly nested difference operators and least squares fits. 

However, the theory is quite involved. The reader is referred to Sjogreen and Yee for the exact 
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formula and the references cited for background. It is noted that the dual purpose wavelet 
sensors (dynamic numerical dissipation controls and grid adaptation indicators) can be a stand 
alone option for a variety of schemes other than what is discussed here. 

It is emphasized here that neither ACM nor wavelet sensors will be able to improve the 
accuracy at the shock and shear locations over the inherent shock-capturing capability of the 
nonlinear dissipation. The accuracy of the shock and shear is dictated by the chosen nonlinear 
dissipation. The role of the sensors is to allow the full amount of numerical dissipation in shock 
and shear regions, and to limit the amount of numerical dissipation in regions immediately 
away from shock and shear locations and the rest of the flow field. Therefore, with a suitable 
sensor, one does not have to use CPU- intensive high order high-resolution shock-capturing 
numerical dissipation, since this type of dissipation generally gives a slightly more accurate 
solution away from discontinuities but exhibits similar shock and shear resolutions as second 
or third-order high-resolution numerical dissipations. 

Full Discretizations: If a multistage time discretization such as the Runge-Kutta method is 
desired, the high order non-dissipative spatial differencing base scheme is applied at every stage 
of the Runge-Kutta method. If viscous terms are present, they use the same order and type of 
base scheme as for the inviscid terms. There are two methods for applying the characteristic 
filter. Method 1 is to apply the filter at every stage of the Runge-Kutta step. Method 2 is 
to apply the filter at the end of the full Runge-Kutta step. For inviscid and strong shock 
interactions, method 1 might be more stable. 

If one desires a time discretization that belongs to the class of linear multistep methods 
(LMMs), e.g.. trapezoidal rule or three-point backward differentiation, then the filter can be 
applied as a numerical dissipation vector in conjunction with the base scheme. The filter in 
this case is evaluated at for explicit LMMs. For implicit LMMs additional similar filters 
evaluated at the n^l time level might be involved. Alternatively, method 2 can be applied to 
LMMs as well. In this case, we apply the filter after the completion of the implicit time step. 

As an example, we illustrate the complete form of the schemes for Runge-Kutta methods 
w'ith the filters applied at the completion of a full Runge-Kutta time step. Let 
solution after one full Runge-Kutta time step using a non-dissipative spatial base scheme. Note 
that if entropy splitting is employed, the base scheme is applied to the split form of the inviscid 
flux derivatives. Then the solution at the next time level is 


rrn-hl 


At 

Ax 


-w 

At 
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are evaluated at 
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III. Numerical Examples 

We summarize the performance of this class of schemes by illustrating four perfect gas test 
cases with distinct flow' properties. The first is inviscid and the last three are compressible full 
Navier-Stokes computations. The four test cases are: (1) a 2-D inviscid horizontally convecting 
vortex with periodic boundary conditions (BCs), (2) a 2-D vortex pairing in a time-developing 
mixing layer with shock waves formed around the vortices, (3) a 2-D shock wave impinging 
on a spatially evolving mixing layer where the evolving vortices must pass through a shock 
wave, which in turn is deformed by the vortex passage, and (4) a direct numerical simulation 
of a 3-D shock-free compressible turbulent plane channel flow. Figures 1-4 show the schematic, 
flow conditions, grid size and the computational domains of the four test cases. 
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In all of the computations the classical fourth-order Runge-Kutta time discretization, and 

the non-compact central spatial interior schemes with the same order of accuracy and type 

of base scheme for the inviscid and viscous terms (if viscosities are present) are employed. If 

filters are used, they are applied at the end of the full Runge-Kutta time step. Roe’s (1981) 

average states are used in (4), along with the Harten and Yee (Yee 1989, Yee et al. 1990) 

second-order upwind TVD dissipation portion for (^)- The parameters p and m are 

2 

set to 1 and a small value of 10“® is added to the denominator of (7) to avoid an extra logical 
statement for the ACM sensor. These will be notated as ACM or WAV (depending on whether 
ACM or wavelet is used as sensor) with the following numbers indicating the order of the 
interior base scheme for the inviscid and viscous terras. For example, ACM66 (WAV66) means 
the use of sixth-order central as the base scheme for both the inviscid and viscous terms, and 
ACM as sensor (wavelet as sensor). In order not to introduce additional notation, inviscid 
flow simulations are designated by the same notation, with the viscous terms not activated. 
Computations using entropy splitting are indicated by adding the letters “ENT” at the end 
as in ACM66-ENT (WAV66-ENT). Computation using 5^. , ^ the sensor is turned off 

and the upwind TVD dissipation portion is used as the filter, will be notated by TVD as in 
TVD66. To examine the performance of the entropy splitting schemes where shock waves are 
absent, the computations also employ only the non-dissipative central base schemes (without 
the filters) designated as CEN66 and CEN66-ENT for sixth order unsplit and split forms of 
the inviscid flux derivatives respectively. The fifth-order WENO scheme of Shu (1997) will be 
notated as WEN05. 

For the Navier-Stokes computations involving entropy splitting, the splitting is applied to 
the inviscid flux terms, and the symmetric form of the viscous flux is not used. Except for the 
periodic BCs and the 3-D channel problem, compatible stable boundary schemes of Strand, and 
Olsson are not used in order to study the benefit of the entropy splitting without considering 
the accompanying stable boundary schemes as a complete package for stability consideration. 
For the 3-D channel flow, a variant of the Strand boundary scheme developed by Carpenter 
et al. (1998) is employed. For this problem, we also reformulate the viscous term to minimize 
odd and even decoupling phenomena before the application of the central scheme (Sandham 
& Yee). For the second and third test cases, w’e lower the order of the base scheme near the 
boundary points for the boundary scheme. For the third test case, for simplicity, slip wall 
BC is used for the lower wall, and the upper y-direction physical BC is overspecified and 
nonreflecting BC is not used. The various ways of imposing physical and boundary schemes 
for the different cases illustrate the robustness of the schemes and are an indication of the 
added value of entropy splitting used in conjunction with interior and boundary schemes that 
do not satisfy the generalized energy estimate. In order to assess the true performance of 
the algorithm, no attempt is made to enhance the resolution using appropriate adaptive grid 
procedures- 

Figures 5-16 show the following comparisons: 

1. Numerical Dissipation Sensor Control vs. No Sensor: Comparisons of TVD66 and ACM66 
for test cases 2 and 3 are shown in Figures 5 and 6. For test case 2, the ACM66 using a 
41 X 41 grid produces a solution slightly more accurate (around the shear regions) than 
TVD66 using a 101 x 101 grid. Note that the widths of the shock resolution of ACM66 and 
TVD66 are the same as discussed in Section II (the paragraph before the full discretization 
subsection). For the TVD66, ACM66 and WAV66 comparison, results should be compared 
with Figures 12 and 13a results. 


10 


2. Entropy Splitting vs. Unsplit: Comparisons of CEN66 vs. CEN66-ENT, and ACM66 
vs. ACM66-ENT for test case 1 are shown in Figures 7-10. The CEN66 diverged shortly 
after f = 55 as opposed to over 6 times the integration time longer for CEN66-ENT 
with an almost perfect vortex preservation. For even longer time integrations, numerical 
dissipation is needed. The amount for ACM66-ENT is less than for ACM66. For all 
simulations involving entropy splitting for test case 1, we set the splitting parameter 
/3 = 1 and p* in (3) to be a constant (for isentropic flow). It seems that ACM66-ENT 
requires a factor of at least less numerical dissipation over the ACM66. This factor is 
precisely the conservative portion of the entropy splitting of the inviscid flux derivatives. 
This factor is, however, greatly reduced if high shear or shock waves are present. 

3. Performance of the 5th-Order WENO Scheme: We compare WEN05 with ACM66 for 
test cases 3 and 4, and various mixing layer problems, including hypersonic speed. In all 
cases. WEN05 is more diffusive than ACM66. Figure 11 shows their comparison for test 
case 3, where the normalized temperature contours (top two) and Mach contours (bottom 
two) at t=113.16 on a 320 x 81 grid for WEN05 and a 321 x 81 grid for ACM66 are 
shown. The vortices are more diffusive in the WEN05 computations. There is a minor 
difference on the two simulations. The WEN05 code has a built-in nonreflecting BC on 
the upper y-direction. It is too costly to fully compare the WEN05 with the 3-D channel 
DNS computations where a fully developed turbulence statistics grid refinement study 
was illustrated in Sandham & Yee (2000) using CEN44-ENT. For WEN05, a smaller 
CFL number is used for stability. We compared the CPU time per time step and accuracy 
in a loose sense with CEN44-ENT. We use the solution at t = 150 from the CEN44-ENT 
simulation as the starting solution for WEN05 and integrate for a time length of 20 and 
compare the turbulence statistics with the CEN44-ENT simulation. Our preliminary study 
shows that WEN05 costs a factor of six more CPU time with less accurate turbulence 
statistics than CEN44-ENT. See Hadjadj et al. (2001) for details. 

4. ACM sensor vs. Wavelet sensor: Sample computations using WAV66 for test cases 2 and 

3 are shown in Figures 12-14. The accuracy of the two wavelet sensors, B-spline wavelet 
(WAV66-BS) or the redundant form of Harten w^avelet (WAV66-RH) for test cases 1-3 
(results not shown) is very similar and the effect on accuracy of the choice of the physical 
vector (density and/or pressure, characteristic variables, or entropy variables W) to 

be sensed is not pronounced. In all cases, no physical problem-dependent parameter has 
to be tuned. The accuracy compared very well with that of the corresponding best tuned 
K for the individual test cases 1-3. In particular, the same accuracy was sustained, using 
WAV66-ENT for either the B-spline or the redundant form of Harten wavelet sensors for 
long time integrations of the vortex convection problems, as ACM66-ENT using k = 0.01 
and — 0.01. Figure 12 shows the test case 2 results which should be compared with 
the ACM66 results in Figure 5. Figures 13 and 14 show the test case 3 results. Figure 13a 
shows the density and pressure contours computations by WAV66-RH using the density 
and pressure as the functions to be sensed (see Figure 6 for comparison with ACM66). 
Figure 13b show’s the corresponding estimated Lipschitz exponent (alpha) for the density 
and pressure at t=T20 and the wavelet sensor itself. Figure 14a shows the wavelet sensor 
applied to the density and pressure in the x and y-directions. and the square root of the 
sum of these quantities in the x and y-directions. Figure 14b shows the corresponding 
contours using the ACM sensor. The wavelet sensor was able to extract the full feature 
of the flow structure far better than the ACM sensor. See Sjogreen and Yee for the exact 
formula and detailed numerical experiments. 
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5. Fourth-Order Central Differencing with Entropy Splitting vs. Spectral Method: A grid 
refinement study of a compressible turbulent DNS and a compcirison of CEN44-ENT 
with the incompressible simulation using a spectral method for test case 4 are shown in 
Figures 15 and 16. For the physical problems, background and grid refinement study, 
see Sandhain k Yee for details. These computations require very long time integrations 
of the full Navier-Stokes equations before fully developed turbulence statistics can be 
obtained. A very 2 u:curate data base is available for the 3-D DNS incompressible turbulence 
simulations. For compcirison, we use a Mach number of 0.1 and a Reynolds number of 180. 
Very accurate fully developed turbulence statistics were obtained using coarse to moderate 
grid sizes and without filters. The results compared well with the best spectral method of 
incompressible Navier-Stokes flow formulations which may use de-aliasing, skew-symmetric 
formulation or energy conservation (e.g., variables on a staggered grid) to obtain a robust 
method. The angle brackets < > denote averages over the homogeneous spatial directions 
and time while double primes denote deviation from mass-weighted (Favre) averages. 

Concluding Remarks 

In addition to the systematic theoretical development of simple parallelizable nonlinear stable 
schemes for IBVPs, their extension and their abilities to accurately simulate a wide range of 
physical applications for the Euler and Navier-Stokes equations, there are two key items that 
might have potential impact on a wide range of complex practical computations. One is the 
possible applicability of entropy splitting beyond the boundary of the theory indicated and the 
other is the wavelet sensors. 

The numerical results indicate that entropy splitting alone can improve nonlinear stability 
even when one employs boundary schemes that do not satisfy the discrete generalized energy 
estimate. This stability property is valuable not just to the class of schemes in question, 
but might also be applicable to other schemes commonly used in practical CFD applications. 
This finding actually is not surprising since one should, if possible, condition the governing 
equations before the application of a suitable numerical method that can handle the type of 
physics inherent in the physical model. 

The proposed wavelet sensors, unlike the ACM sensor, can detect most of the distinct flow 
features, including turbulence, leading to an automatic selection of the appropriate distribution 
of numerical dissipation. In addition, these wavelet sensors are free of physical problem- 
dependent parameters for the three test cases, and they are also good grid adaptation indicators 
(Gerritsen k Olsson) when compared to the ones commonly used in practice. Consequently, a 
new dual purpose adaptive method is readily available leading to dynamic numerical dissipation 
controls and improved grid adaptation indicators. This dual purpose adaptive method can 
also serve as a stand alone option for other numerical schemes. These are subjects of ongoing 
research. 
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Tables Captions 

Table 1. Flow Chart of the efficient low dissipative high order schemes. 

Figures Captions 

Fig- 1 

Fig. 2 
Fig. 3 

Fig. 4 

Fig. 5. Vortex pairing: Comparison of ACM66 with TVD66 = 1) with reference 

solution ACM44 (201 x 201 grid), illustrated by the normalized temperature contours 
at time t = 160 on 101 x 101 and 41 x 41 grids with n = 0.7 for the nonlinear fields 
and K = 0.35 for the linear fields. 

Fig. 6. Shock-shear layer interaction: Comparison of ACM66 with TVD66 = 1) with 

reference solution ACM44 (641 x 161 grid), illustrated by the density contours at 
t = 120 on a 321 X 81 grid with k, = 0.35 for the nonlinear fields and k = 0.175 for the 
linear fields. 

Fig. 7. Convecting Vortex: Comparison of CEN66 with the exact solution (I.C.), illustrated 
by density contours t = 20, 30, 40, 50, 55 for At = 0.01 on a 80 x 79 grid. 

Fig. 8. Convecting Vortex: Comparison of CEN66-ENT with the exact solution (I.C.). illus- 
trated by density contours at t — 100,200,300,400,500 for At = 0.01 on a 80 x 79 
grid. 

Fig. 9. Convecting Vortex: Comparison of ACM66 with the exact solution (I.C.), illustrated 
by density contours at t = 100,200,300,400,500, 600,700,800,900,1000,1100 for 
At = 0.04 and k, = 0.04 on a 80 x 79 grid. 

Fig. 10. Convecting Vortex: Comparison of ACM66-ENT with the exact solution (I.C.), illus- 
trated by density contours at t = 100,200,300,400,500, 600,700,800,900,1000,1100 
for At — 0.01 and k = 0.01 on a 80 x 79 grid. 

Fig. 11. Shock-shear layer interaction: Comparison of ACM66 with WEN05, illustrated by 
the normalized temperature contours (top two plots) and Mach contours (bottom two 
plots) at t — 113.16 on a 321 x 81 grid for ACM66 with k = 0.35 for the nonlinear 
fields and ac = 0.175 for the linear fields, and a 320 x 81 for the WEN05. 

Fig. 12. Vortex pairing: WAV66-RH computations, illustrated by the density and normalized 
temperature contours at time t = 160 on a 101 x 101 grid using the density and 
pressure as functions to be sensed. 

Fig. 13. Shock-shear layer interaction: WAV66-RH computations illustrated by (a) density 
and pressure contours (top 2 plots), and (b) estimated Lipschitz exponent a and the 
wavelet sensor itself ibottom 2 plots) using the density and pressure as functions to 
be sensed at ^ = 120 on a 321 x 81 grid. 

Fig. 14. Shock-shear layer interaction: Comparison of ACM66 with WAV66 at ^ = 120 on a 
321 x 31 grid with k = 0.35 for the nonlinear fields and k = 0.175 for the linear fields 
for ACM66, illustrated by (a) contours of the sensor used by WAV66-RH applied to 
the density and pressure in the x and y-directions, and the square root of the sum 
of these quantity (top 3 plots), and (b) the corresponding contours using the ACM 
sensor (bottom 3 plots). 

Fig. 15. 3-D channel: Effect of grid refinement of a compressible turbulent DNS on (a) root 
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mean square of u, v and w and (b) Reynolds stress and total stress, using CEN44-ENT 
atnd Carpenter et al. boundary scheme. 

Fig. 16. 3-D channel: Comparison of root mean square of u, u, and it; on a 30 x 30 x 101 grid 
of a compressible turbulent DNS using CEN44-ENT and Carpenter et al. boundciry 
scheme with an incompressible turbulent DNS result on a 32 x 32 x 81 grid using 
spectral method. 
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Table 1. Flow Chart of the efficient low dissipative high order schemes. 
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Fig. 2 
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Fig. 3 
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Fig. 5. Vortex pairing: Comparison of ACM66 with TVD66 («0' , = 1) with reference 

solution ACM44 (201 x 201 grid), illustrated by the normalized temperature contours 
at time t = 160 on 101 x 101 and 41 x 41 grids with k = 0.7 for the nonlinear fields 
and K = 0.35 for the linear fields. 


Shock Impingement on Mixing Layer 

Density Contours at t=120 



Fig. 6. Shock-shear layer interaction: Comparison of ACM66 with TVD66 = 1 ) with 

reference solution ACM44 (641 x 161 grid), illustrated by the density contours at 
t = 120 on a 321 x 81 grid with « = 0.35 for the nonlinear fields and k = 0.175 for the 
linear fields. 
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Fig. 8. Convecting Vortex: Comparison of CEN66-ENT with the exact solution (I.C.), illus- 
trated by density contours at t = 100,200,300,400,500 for At = 0.01 on a 80 x 79 
grid. 













Fig. 9. Convecting Vortex: Compajison of ACM66 with the exact solution (I.C.), illustrated 
by density contours at t = 100,200,300,400,500, 600,700,800,900,1000,1100 for 
At = 0.04 and « = 0.04 on a 80 x 79 grid. 
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Fig. 10. Convecting Vortex; Comparison of ACM66-ENT with the exact solution (I.C.), illus- 
trated by density contours at t = 100,200,300,400,500, 600,700,800,900,1000,1100 
for A< = 0.01 and k = 0.01 on a 80 x 79 grid. 
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ACM66 vs. WEN05, 321 x 81 grid, Ta113.16 



Fig. 11. Shock-shear layer interaction: Comparison of ACM66 with WEN05, illustrated by 
the normalized temperature contours (top two plots) and Mach contours (bottom two 
plots) at t = 113.16 on a 321 x 81 grid for ACM66 with k = 0.35 for the nonlinear 
,1 ■ n 17T, f,,,, ,r,„ floPk 'i>vl ‘^‘^n v «l fr„- 




Vortex Pairing at Mg = 0.8 

WAV66-RH, t=160, 101 x 101 grid 



Fig. 12. Vortex pairing; WAV66-RH computations, illustrated by the density and normalized 
temperature contours at time t = 160 on a 101 x 101 grid using the density and 
pressure as functions to be sensed. 
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Fig. 13. Shock-shecir layer interaction: WAV66-RH computations illustrated by (a) density 
and pressure contours, and (b) estimated Lipschitz exponent q and the wavelet sensor 
itself using the density and pressure as functions to be sensed at t = 120 on a 321 x 81 
grid. 
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WAV66-RH, t=120, 321 x 81 grid 
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Fig. 14. Shock-shear layer interaction; Comparison of ACM66 with WAV66 at t = 120 on a 
321 X 81 grid with k = 0.35 for the nonlinear fields and k = 0.175 for the linear fields 
for ACM66, illustrated by (a) contours of the sensor used by WAV66-RH applied to 
the density and pressure in the x and y-directions, and the square root of the sum of 
these quantity, and (b) the corresponding contours using the ACM sensor. 
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Fig. 15. 3-D channel: Effect of grid refinement of a compressible turbulent DNS on (a) root 
mean square of u, v and w and (b) Reynolds stress and total stress, using CEN44-ENT 
and Carpenter et al. boundary scheme. 



Fig. 16. 3-D channel: Comparison of root mean square of u, u, and tx; on a 30 x 30 x 101 grid 
of a compressible turbulent DNS using CEN44-ENT and Carpenter et al. boundary 
scheme with an incompressible turbulent DNS result on a 32 x 32 x 81 grid using 
spectral method. 






